      subroutine sanmjin

      use param

#include <global.h>
      include 'ntbytes.h'
      include 'geom.h'
      include 'xsec.h'
      include 'sanm.h'
      include 'sanmjin.h'
      include 'eigv.h'
      include 'setls.h'
      include 'cmfd.h'

      real(NBF) cp011,cp012,cp021,cp022,cp111,cp112,cp121,cp122,cp211,cp212 &
     ,         cp221,cp222,cp311,cp312,cp321,cp322,cp411,cp412,cp421,cp422 &
     ,         krt11,krt12,krt21,krt22,rkrt11,rkrt12,rkrt21,rkrt22 &
     ,         ksq11,ksq12,ksq21,ksq22,rksq11,rksq12,rksq21,rksq22 &
     ,         skrt11,skrt12,skrt21,skrt22,ckrt11,ckrt12,ckrt21,ckrt22 &
     ,         eig11,eig12,eig21,eig22,beta11,beta12,beta21,beta22 &
     ,         rskrt11,rskrt12,rskrt21,rskrt22,j11_err,j12_err,j21_err,j22_err &
     ,         j11_old,j12_old,j21_old,j22_old,jpl,jpr,ccpl11,ccpl12,ccpl21,ccpl22 &
     ,         dccpl11,dccpl12,dccpl21,dccpl22,ccpr11,ccpr12,ccpr21,ccpr22 &
     ,         dccpr11,dccpr12,dccpr21,dccpr22,tau11,tau12,tau21,tau22,tau31,tau32 &
     ,         tau41,tau42,gamma1,gamma2,mu11,mu12,mu21,mu22,mu31,mu32,mu41,mu42 &
     ,         mu51,mu52,zeta11,zeta12,zeta21,zeta22,zeta31,zeta32,zeta41,zeta42 &
     ,         zeta51,zeta52,zeta61,zeta62,zeta71,zeta72,zeta81,zeta82,yeta11,yeta12 &
     ,         yeta21,yeta22,yeta31,yeta32,yeta41,yeta42,yeta51,yeta52,yeta61,yeta62 &
     ,         lm11,lm12,lm21,lm22,lm31,lm32,lc1,lc2,lc3,rm11,rm13,rm21,rm23,rm31,rm33 &
     ,         rc1,rc2,rc3,m11,m12,m13,m21,m22,m23,m31,m32,m33,bb1,bb2,bb3
      real(NBF),dimension(ng) :: af21,af22

! Update the source coefficients with the average source from CMFD
#define srcup
#ifdef srcup
      do im=1,ng
        do i=1,tnode
          do k=1,4
            if (qs(0,i,im).ne.0) then
              qs(k,i,im)=src(1,i,im)*(qs(k,i,im)/qs(0,i,im))
            endif              
          enddo        
          qs(0,i,im)=src(1,i,im)
        enddo
      enddo
#endif

!define test
#ifdef test
      qs(0,1,1)=0.299987
      qs(1,1,1)=0
      qs(2,1,1)=-0.0608546
      qs(3,1,1)=0
      qs(4,1,1)=-0.0237439
      qs(0,2,1)=0.270686
      qs(1,2,1)=-0.109798
      qs(2,2,1)=0.0659352
      qs(3,2,1)=0.0027707
      qs(4,2,1)=0.0609383
      avgflx(1,1,1)=10.38437
      avgflx(1,2,1)=7.718542
      avgflx(2,1,1)=-0.00706014
      avgflx(2,2,1)=-0.0375938
      jin(1,1,1)=2.6208631129358784
      jin(2,2,1)=0.44787943904118926

      qs(0,1,2)=0.21215
      qs(1,1,2)=0
      qs(2,1,2)=-0.00422986
      qs(3,1,2)=0
      qs(4,1,2)=0.00273227
      qs(0,2,2)=0.122439
      qs(1,2,2)=-0.0556648
      qs(2,2,2)=-0.00992264
      qs(3,2,2)=0.000446542
      qs(4,2,2)=-0.00197414
      avgflx(1,1,2)=2.1103759
      avgflx(1,2,2)=0.5986204
      avgflx(2,1,2)=-0.0058428
      avgflx(2,2,2)=0.0057008
      jin(1,1,2)=0.241190
      jin(2,2,2)=0.303246

      eigv=0.956867
      reigv=1/eigv
#endif

! Decide the coefficients of the homogeneous solution          
!      do k=1,2
!        do i=k,tnode-1,2
         do i=1,tnode-1   
          j=i+1
          j11_old=0
          j12_old=0
          j21_old=0
          j22_old=0
          do iin=1,1
            fsc=0
            do im=1,ng
! substitutions    
              krt11=krt(1,i,im)
              krt12=krt(1,j,im)
              krt21=krt(2,i,im)
              krt22=krt(2,j,im)
              ksq11=ksq(1,i,im)
              ksq12=ksq(1,j,im)
              ksq21=ksq(2,i,im)
              ksq22=ksq(2,j,im)
              rkrt11=rkrt(1,i,im)
              rkrt12=rkrt(1,j,im)
              rkrt21=rkrt(2,i,im)
              rkrt22=rkrt(2,j,im)
              rksq11=rksq(1,i,im)
              rksq12=rksq(1,j,im)
              rksq21=rksq(2,i,im)
              rksq22=rksq(2,j,im)

              skrt11=sinh(krt11)
              skrt12=sinh(krt12)
              skrt21=sinh(krt21)
              skrt22=sinh(krt22)
              ckrt11=cosh(krt11)
              ckrt12=cosh(krt12)
              ckrt21=cosh(krt21)
              ckrt22=cosh(krt22)

              eig11=eigvec1(i,im)
              eig12=eigvec1(j,im)
              eig21=eigvec2(i,im)
              eig22=eigvec2(j,im)
              beta11=betax(1,i,im)
              beta12=betax(1,j,im)
              beta21=betax(2,i,im)
              beta22=betax(2,j,im)

              rskrt11=1/sinh(krt11)
              rskrt12=1/sinh(krt12)
              rskrt21=1/sinh(krt21)
              rskrt22=1/sinh(krt22)
! Incoming partial current
              jpl=jin(1,i,im)
              jpr=jin(2,j,im)
! Decide the coefficients of the particular solution
! cp/order/moment/l-rnode/
              cp011=imq1(i,im)*(5*qs(4,i,im)*(21+2*ksq11)+qs(0,i,im)*kqu(1,i,im) &
                   +3*qs(2,i,im)*ksq11)*rksx(1,i,im)
              cp012=imq1(j,im)*(5*qs(4,j,im)*(21+2*ksq12)+qs(0,j,im)*kqu(1,j,im) &
                   +3*qs(2,j,im)*ksq12)*rksx(1,j,im)
              cp021=imq2(i,im)*(5*qs(4,i,im)*(21+2*ksq21)+qs(0,i,im)*kqu(2,i,im) &
                   +3*qs(2,i,im)*ksq21)*rksx(2,i,im)
              cp022=imq2(j,im)*(5*qs(4,j,im)*(21+2*ksq22)+qs(0,j,im)*kqu(2,j,im) &
                   +3*qs(2,j,im)*ksq22)*rksx(2,j,im)
              cp111=imq1(i,im)*(qs(1,i,im)*ksq11+15*qs(3,i,im))*rkqu(1,i,im)
              cp112=imq1(j,im)*(qs(1,j,im)*ksq12+15*qs(3,j,im))*rkqu(1,j,im)
              cp121=imq2(i,im)*(qs(1,i,im)*ksq21+15*qs(3,i,im))*rkqu(2,i,im)
              cp122=imq2(j,im)*(qs(1,j,im)*ksq22+15*qs(3,j,im))*rkqu(2,j,im)
              cp211=imq1(i,im)*(qs(2,i,im)*ksq11+35*qs(4,i,im))*rkqu(1,i,im)
              cp212=imq1(j,im)*(qs(2,j,im)*ksq12+35*qs(4,j,im))*rkqu(1,j,im)
              cp221=imq2(i,im)*(qs(2,i,im)*ksq21+35*qs(4,i,im))*rkqu(2,i,im)
              cp222=imq2(j,im)*(qs(2,j,im)*ksq22+35*qs(4,j,im))*rkqu(2,j,im)
              cp311=imq1(i,im)*qs(3,i,im)*rksq11
              cp312=imq1(j,im)*qs(3,j,im)*rksq12
              cp321=imq2(i,im)*qs(3,i,im)*rksq21
              cp322=imq2(j,im)*qs(3,j,im)*rksq22
              cp411=imq1(i,im)*qs(4,i,im)*rksq11
              cp412=imq1(j,im)*qs(4,j,im)*rksq12
              cp421=imq2(i,im)*qs(4,i,im)*rksq21
              cp422=imq2(j,im)*qs(4,j,im)*rksq22
! left node
              ccpl11=cp011-cp111+cp211-cp311+cp411
              ccpl21=cp021-cp121+cp221-cp321+cp421
              dccpl11=cp111-3*cp211+6*cp311-10*cp411 
              dccpl21=cp121-3*cp221+6*cp321-10*cp421 
              ccpl12=cp011+cp111+cp211+cp311+cp411
              ccpl22=cp021+cp121+cp221+cp321+cp421
              dccpl12=cp111+3*cp211+6*cp311+10*cp411 
              dccpl22=cp121+3*cp221+6*cp321+10*cp421 
! right node
              ccpr11=cp012-cp112+cp212-cp312+cp412
              ccpr21=cp022-cp122+cp222-cp322+cp422
              dccpr11=cp112-3*cp212+6*cp312-10*cp412
              dccpr21=cp122-3*cp222+6*cp322-10*cp422
              ccpr12=cp012+cp112+cp212+cp312+cp412
              ccpr22=cp022+cp122+cp222+cp322+cp422
              dccpr12=cp112+3*cp212+6*cp312+10*cp412
              dccpr22=cp122+3*cp222+6*cp322+10*cp422
!
              tau11=beta11*(2+eig11)
              tau12=beta12*(2+eig12)
              tau21=beta11*(2+eig21)
              tau22=beta12*(2+eig22)
              tau31=5+4*eig11
              tau32=5+4*eig12
              tau41=5+4*eig21
              tau42=5+4*eig22
              gamma1=1/(eig11-eig21)
              gamma2=1/(eig12-eig22)
              mu11=-ckrt11*krt11*tau11
              mu12=-ckrt12*krt12*tau12
              mu21=-ckrt21*krt21*tau21
              mu22=-ckrt22*krt22*tau22
              mu31=-0.0625*skrt11*(5+4*eig11)-0.5*ckrt11*krt11*tau11
              mu32=-0.0625*skrt12*tau32-0.5*ckrt12*krt12*tau12
              mu41=-0.0625*skrt21*(5+4*eig21)-0.5*ckrt21*krt21*tau21
              mu42=-0.0625*skrt22*tau42-0.5*ckrt22*krt22*tau22
              mu51=1/(mu11*mu41-mu21*mu31)
              mu52=1/(mu12*mu42-mu22*mu32)
!
              zeta11=-eig21*krt11*rskrt11*gamma1
              zeta12=-eig22*krt12*rskrt12*gamma2
              zeta21=eig11*krt21*rskrt21*gamma1
              zeta22=eig12*krt22*rskrt22*gamma2
              zeta31=-(skrt11*zeta11*krt11*tau11+skrt21*zeta21*krt21*tau21)
              zeta32=-(skrt12*zeta12*krt12*tau12+skrt22*zeta22*krt22*tau22)
              zeta41=-0.5*zeta31+0.0625*(ckrt11*zeta11*tau31+ckrt21*zeta21*tau41)
              zeta42=-0.5*zeta32+0.0625*(ckrt12*zeta12*tau32+ckrt22*zeta22*tau42)
              zeta51=mu51*(mu21*zeta41-mu41*zeta31)
              zeta52=mu52*(mu22*zeta42-mu42*zeta32) !
              zeta61=mu51*(mu31*zeta31-mu11*zeta41)
              zeta62=mu52*(mu32*zeta32-mu12*zeta42) !
              zeta71=mu41*mu51
              zeta72=mu42*mu52
              zeta81=-mu31*mu51
              zeta82=-mu32*mu52
!
              yeta11=krt11*rskrt11*(avgflx(1,i,im)*gamma1-cp011)
              yeta12=krt12*rskrt12*(avgflx(1,j,im)*gamma2-cp012)
              yeta21=-krt21*rskrt21*(avgflx(1,i,im)*gamma1+cp021)
              yeta22=-krt22*rskrt22*(avgflx(1,j,im)*gamma2+cp022)
              yeta31=-(skrt11*yeta11*krt11*tau11+skrt21*yeta21*krt21*tau21+tau11*dccpl12+tau21*dccpl22)
              yeta32=skrt12*yeta12*krt12*tau12+skrt22*yeta22*krt22*tau22-tau12*dccpr11-tau22*dccpr21
              yeta41=0.5*(skrt11*yeta11*krt11*tau11+skrt21*yeta21*krt21*tau21-tau11*dccpl11-tau21*dccpl21) &
                   +0.0625*(ckrt11*yeta11*tau31+ckrt21*yeta21*tau41+tau31*ccpl11+tau41*ccpl21)
              yeta42=0.5*(skrt12*yeta12*krt12*tau12+skrt22*yeta22*krt22*tau22+tau12*dccpr12+tau22*dccpr22) &
                   +0.0625*(ckrt12*yeta12*tau32+ckrt22*yeta22*tau42+tau32*ccpr12+tau42*ccpr22)
              yeta51=mu51*(mu21*yeta41-mu21*jpl-mu41*yeta31)
              yeta52=mu52*(mu22*jpr-mu42*yeta32-mu22*yeta42) !
              yeta61=mu51*(mu11*jpl+mu31*yeta31-mu11*yeta41)
              yeta62=mu52*(mu12*yeta42+mu32*yeta32-mu12*jpr) !
! Coefficients of the Cosh term
              cosa(1,i,im)=zeta11*avgflx(2,i,im)+yeta11    
              cosa(1,j,im)=zeta12*avgflx(2,j,im)+yeta12         
              cosa(2,i,im)=zeta21*avgflx(2,i,im)+yeta21
              cosa(2,j,im)=zeta22*avgflx(2,j,im)+yeta22
! coeff. matrix of J, flux2 and constant/ left node (3X3)
              lm11=-beta21*(ckrt11*zeta71*krt11+ckrt21*zeta81*krt21)
              lm12=-beta21*(krt11*(skrt11*zeta11+ckrt11*zeta51)+krt21*(skrt21*zeta21+ckrt21*zeta61))
              lm21=skrt11*eig11*zeta71+skrt21*eig21*zeta81
              lm22=eig11*(ckrt11*zeta11+skrt11*zeta51)+eig21*(ckrt21*zeta21+skrt21*zeta61)
              lm31=skrt11*zeta71+skrt21*zeta81
              lm32=ckrt11*zeta11+ckrt21*zeta21+skrt11*zeta51+skrt21*zeta61
              lc1=-beta21*(skrt11*yeta11*krt11+skrt21*yeta21*krt21+ckrt11*krt11*yeta51 &
                 +ckrt21*krt21*yeta61+dccpl12+dccpl22)
              lc2=eig11*(ckrt11*yeta11+ccpl12+skrt11*yeta51)+eig21*(ckrt21*yeta21+ccpl22+skrt21*yeta61)
              lc3=ckrt11*yeta11+ckrt21*yeta21+skrt11*yeta51+skrt21*yeta61+ccpl12+ccpl22
! right node (3X3)
              rm11=-beta22*(ckrt12*zeta72*krt12+ckrt22*zeta82*krt22)
              rm13=beta22*(krt12*(skrt12*zeta12+ckrt12*zeta52)+krt22*(skrt22*zeta22+ckrt22*zeta62)) 
              rm21=-(skrt12*eig12*zeta72+skrt22*eig22*zeta82)
              rm23=eig12*(ckrt12*zeta12+skrt12*zeta52)+eig22*(ckrt22*zeta22+skrt22*zeta62) 
              rm31=-(skrt12*zeta72+skrt22*zeta82)
              rm33=ckrt12*zeta12+ckrt22*zeta22+skrt12*zeta52+skrt22*zeta62 
              rc1=beta22*(skrt12*yeta12*krt12+skrt22*yeta22*krt22-ckrt12*krt12*yeta52 &
                 -ckrt22*krt22*yeta62-dccpr11-dccpr21)
              rc2=eig12*(ckrt12*yeta12-skrt12*yeta52+ccpr11)+eig22*(ckrt22*yeta22-skrt22*yeta62+ccpr21)
              rc3=ckrt12*yeta12+ckrt22*yeta22-skrt12*yeta52-skrt22*yeta62+ccpr11+ccpr21
!
              m11=lm11-rm11
              m12=lm12
              m13=-rm13
              m21=lm21-rm21
              m22=lm22
              m23=-rm23
              m31=lm31-rm31
              m32=lm32
              m33=-rm33
              bb1=rc1-lc1
              bb2=rc2-lc2
              bb3=rc3-lc3  
! net current J
              jnet(1,j,im)=(bb3*(-m13*m22+m12*m23)+bb2*(m13*m32-m12*m33)+bb1*(-m23*m32+m22*m33)) &
                         /(m13*(-m22*m31+m21*m32)+m12*(m23*m31-m21*m33)+m11*(-m23*m32+m22*m33))
! 2nd moment flux
              avgflx(2,i,im)=(bb3*(m13*m21-m11*m23)+bb2*(-m13*m31+m11*m33)+bb1*(m23*m31-m21*m33)) &
                            /(m13*(-m22*m31+m21*m32)+m12*(m23*m31-m21*m33)+m11*(-m23*m32+m22*m33))
              avgflx(2,j,im)=(bb3*(-m12*m21+m11*m22)+bb2*(m12*m31-m11*m32)+bb1*(-m22*m31+m21*m32)) &
                            /(m13*(-m22*m31+m21*m32)+m12*(m23*m31-m21*m33)+m11*(-m23*m32+m22*m33)) 
! 2nd moment current J2
              jnet(2,j,im)=rm11*jnet(1,j,im)+rm13*avgflx(2,j,im)+rc1
! Coefficients of the Sinh term
              sinb(1,i,im)=zeta51*avgflx(2,i,im)+zeta71*jnet(1,j,im)+yeta51
              sinb(1,j,im)=-zeta52*avgflx(2,j,im)+zeta72*jnet(1,j,im)+yeta52
              sinb(2,i,im)=zeta61*avgflx(2,i,im)+zeta81*jnet(1,j,im)+yeta61
              sinb(2,j,im)=-zeta62*avgflx(2,j,im)+zeta82*jnet(1,j,im)+yeta62
! Surface flux
              if (i.eq.1) then
                sflux(1,i,im)=eig11*(ckrt11*cosa(1,i,im)-skrt11*sinb(1,i,im)+ccpl11) &
                             +eig21*(ckrt21*cosa(2,1,im)-skrt21*sinb(2,i,im)+ccpl21)
                sflux(2,i,im)=ckrt11*cosa(1,i,im)+ckrt21*cosa(2,i,im)-skrt11*sinb(1,i,im) &
                             -skrt21*sinb(2,i,im)+ccpl11+ccpl21
              elseif (j.eq.tnode) then
                sflux(1,j+1,im)=eig12*(ckrt12*cosa(1,j,im)+skrt12*sinb(1,j,im)+ccpr12) &
                               +eig22*(ckrt22*cosa(2,j,im)+skrt22*sinb(2,j,im)+ccpr22)
                sflux(2,j+1,im)=ckrt12*cosa(1,j,im)+ckrt22*cosa(2,j,im)+skrt12*sinb(1,j,im) &
                               +skrt22*sinb(2,j,im)+ccpr12+ccpr22
              endif
              sflux(1,j,im)=eig12*(-skrt12*sinb(1,j,im)+ckrt12*cosa(1,j,im)+ccpr11) &
                           +eig22*(-skrt22*sinb(2,j,im)+ckrt22*cosa(2,j,im)+ccpr21)
              sflux(2,j,im)=ckrt12*cosa(1,j,im)+ckrt22*cosa(2,j,im) &
                           -(skrt12*sinb(1,j,im)+skrt22*sinb(2,j,im))+ccpr11+ccpr21
! incoming current
!              jin(1,j,im)=0.25*sflux(1,j,im)+0.5*jnet(1,j,im)+0.3125*sflux(2,j,im)
!              jin(2,j,im)=0.25*sflux(1,j,im)-0.5*jnet(1,j,im)+0.3125*sflux(2,j,im)
#define srcup2
#ifdef srcup2
! 5 coefficients of the 4-th order flux
              cf(0,i,im)=cosa(1,i,im)*eig11*skrt11*rkrt11+cosa(2,i,im)*eig21*skrt21*rkrt21 &
                      +eig11*cp011+eig21*cp021 
            cf(0,j,im)=cosa(1,j,im)*eig12*skrt12*rkrt12+cosa(2,j,im)*eig22*skrt22*rkrt22 &
                      +eig12*cp012+eig22*cp022
            cf(1,i,im)=3*(sinb(1,i,im)*eig11*(ckrt11*rkrt11-skrt11*rksq11) &
                      +sinb(2,i,im)*eig21*(ckrt21*rkrt21-skrt21*rksq21)) &
                      +eig11*cp111+eig21*cp121
            cf(1,j,im)=3*(sinb(1,j,im)*eig12*(ckrt12*rkrt12-skrt12*rksq12) &
                      +sinb(2,j,im)*eig22*(ckrt22*rkrt22-skrt22*rksq22)) &
                      +eig12*cp112+eig22*cp122
            cf(2,i,im)=5*(cosa(1,i,im)*eig11*rkrt11*(-3*ckrt11*rkrt11+3*skrt11*rksq11+skrt11) &
                      +cosa(2,i,im)*eig21*rkrt21*(-3*ckrt21*rkrt21+3*skrt21*rksq21+skrt21)) &
                      +eig11*cp211+eig21*cp221
            cf(2,j,im)=5*(cosa(1,j,im)*eig12*rkrt12*(-3*ckrt12*rkrt12+3*skrt12*rksq12+skrt12) &
                      +cosa(2,j,im)*eig22*rkrt22*(-3*ckrt22*rkrt22+3*skrt22*rksq22+skrt22)) &
                      +eig12*cp212+eig22*cp222
            cf(3,i,im)=7*(sinb(1,i,im)*eig11*rksq11*(ckrt11*rkrt11*(15+ksq11)-15*skrt11*rksq11-6*skrt11) &
                      +sinb(2,i,im)*eig21*rksq21*(ckrt21*rkrt21*(15+ksq21)-15*skrt21*rksq21-6*skrt21)) &
                      +eig11*cp311+eig21*cp321
            cf(3,j,im)=7*(sinb(1,j,im)*eig12*rksq12*(ckrt12*rkrt12*(15+ksq12)-15*skrt12*rksq12-6*skrt12) &
                      +sinb(2,j,im)*eig22*rksq22*(ckrt22*rkrt22*(15+ksq22)-15*skrt22*rksq22-6*skrt22)) &
                      +eig12*cp312+eig22*cp322
            cf(4,i,im)=9*(cosa(1,i,im)*eig11*rkrt11*(rksq11*(-5*ckrt11*rkrt11*(21+2*ksq11)+105*skrt11*rksq11 &
                      +45*skrt11)+skrt11)+cosa(2,i,im)*eig21*rkrt21*(rksq21*(-5*ckrt21*rkrt21*(21+2*ksq21) &
                      +105*skrt21*rksq21+45*skrt21)+skrt21))+eig11*cp411+eig21*cp421
            cf(4,j,im)=9*(cosa(1,j,im)*eig12*rkrt12*(rksq12*(-5*ckrt12*rkrt12*(21+2*ksq12)+105*skrt12*rksq12 &
                      +45*skrt12)+skrt12)+cosa(2,j,im)*eig22*rkrt22*(rksq22*(-5*ckrt22*rkrt22*(21+2*ksq22) &
                      +105*skrt22*rksq22+45*skrt22)+skrt22))+eig12*cp412+eig22*cp422
!            apsi11=cp011+cosa(1,i,im)*skrt11*rkrt11
!            apsi12=cp012+cosa(1,j,im)*skrt12*rkrt12
!            apsi21=cp021+cosa(2,i,im)*skrt21*rkrt21
!            apsi22=cp022+cosa(2,j,im)*skrt22*rkrt22
!            temp11=cp011*eig11+cp021*eig21+cosa(1,i,im)*eig11*skrt11*rkrt11 &
!                   +cosa(2,i,im)*eig21*skrt21*rkrt21
!            temp12=cp012*eig12+cp022*eig22+cosa(1,j,im)*eig12*skrt12*rkrt12 &
!                   +cosa(2,j,im)*eig22*skrt22*rkrt22
!              temp21=cp011+cp021+cosa(1,i,im)*skrt11*rkrt11+cosa(2,i,im)*skrt21*rkrt21
!              temp22=cp012+cp022+cosa(1,j,im)*skrt12*rkrt12+cosa(2,j,im)*skrt22*rkrt22
! Update the 4-th order source coefficients
              do l=0,4
                fsc(l,i)=fsc(l,i)+xsnf(i,im)*cf(l,i,im)
                fsc(l,j)=fsc(l,j)+xsnf(j,im)*cf(l,j,im)
              enddo
            enddo  ! for im
            do im=1,ng
              scs=0  
              do m=1,ng   
                do l=0,4
                  scs(l,i)=scs(l,i)+xssm(i,im)%from(m)*cf(l,i,m)      
                  scs(l,j)=scs(l,j)+xssm(j,im)%from(m)*cf(l,j,m)      
                enddo
              enddo  
              do l=0,4
                qs(l,i,im)=reigv*xchi(i,im)*fsc(l,i)+scs(l,i)          
                qs(l,j,im)=reigv*xchi(j,im)*fsc(l,j)+scs(l,j)             
              enddo
            enddo ! for im
!            print*,iin
!            print*,jnet(1,j,1),jnet(1,j,2)
!            print*,jnet(2,j,1),jnet(2,j,2)
#endif
! Current Convergence condition
            j11_err=abs(jnet(1,j,1)-j11_old)
            j12_err=abs(jnet(1,j,2)-j12_old)
            j21_err=abs(jnet(2,j,1)-j21_old)
            j22_err=abs(jnet(2,j,2)-j22_old)
            j11_old=jnet(1,j,1)
            j12_old=jnet(1,j,2)
            j21_old=jnet(2,j,1)
            j22_old=jnet(2,j,2)
            if (j21_err.lt.epsm3.and.j22_err.lt.epsm3) then
              if (j11_err.lt.epsm3.and.j12_err.lt.epsm3) then
               print*,"node=",j,"nodal=",iin
               exit
              endif
            endif
          enddo ! for iin
        enddo ! for i
!      enddo ! for k

! surface flux ratio : To prevent a lagged update of the 2nd moment surface flux in the CMFD
      do im=1,ng
        do j=1,tnode+1
!          sfrt(j,im)=sflux(2,j,im)/sflux(1,j,im)
          sfrt(j,im)=sflux(2,j,im)
        enddo
        avgflx(2,0,im)=avgflx(2,1,im)
        avgflx(2,tnode+1,im)=avgflx(2,tnode,im)
      enddo

      end subroutine
